
R version 3.4.0 (2017-04-21) -- "You Stupid Darkness"
Copyright (C) 2017 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

Microsoft R Open 3.4.0
The enhanced R distribution from Microsoft
Microsoft packages Copyright (C) 2017 Microsoft Corporation

Using the Intel MKL for parallel mathematical computing(using 2 cores).

Default CRAN mirror snapshot taken on 2017-05-01.
See: https://mran.microsoft.com/.

[Previously saved workspace restored]

> rm(list=ls())
> set.seed(123456)
> 
> png(filename = "05testplot.png", width=6, height=6.5, units="in", res=400)
> 
> crit <- qt(0.975, df=100)
> 
> x <- seq(from=-1, to=8, by=0.01)
> null.plot <- dt(x, df=100)
> plot(null.plot ~ x, type="n", ylab="probability density", xlab="estimated t statistic",
+      ylim=c(0, 0.55))
> 
> x.shade <- seq(from=crit, to=8, by=0.01)
> y.shade <- dt(x.shade, df=100, ncp=3)
> polygon(c(min(x.shade),x.shade), c(min(y.shade), y.shade), col=gray(0.8, alpha=1), border=NA)
> 
> y.shade <- dt(x.shade, df=100)
> polygon(c(min(x.shade),x.shade), c(min(y.shade), y.shade), col=gray(0.5, alpha=1), border=NA)
> 
> null.plot <- dt(x, df=100)
> lines(null.plot ~ x)
> 
> alt.plot <- dt(x, df=100, ncp=3)
> lines(alt.plot ~ x, lwd=2)
> 
> abline(v = crit, lty=2)
> 
> arrows(5.5, 0.3, 3, 0.15, length=0.1)
> power.value <- round( pt(crit, df=100, ncp=3, lower.tail = F) + pt(-crit, df=100, ncp=3, lower.tail = T), 4)
> text(6, 0.325, labels=(paste( c("power = ",  power.value), collapse = " "  ) ) )
> 
> arrows(5.25, 0.15, 2.5, 0.025, length=0.1)
> size.value <- round( pt(crit, df=100, lower.tail = F) + pt(-crit, df=100, lower.tail = T), 4)
> text(6.5, 0.175, labels=(paste( c("size = ", size.value), collapse = " "  ) ) )
> 
> legend("topright", lty=c(1,1,2), lwd=c(1,2,1), pt.cex=1, cex=0.8,
+        legend=c("t density under null", "t density, ncp = 3", "critical-t, alpha = 0.05"))
> 
> dev.off()
null device 
          1 
> 
> 
> 
> rm(list=ls())
> set.seed(123456)
> 
> png(filename = "005testplot.png", width=6, height=6.5, units="in", res=400)
> 
> crit <- qt(0.9975, df=100)
> 
> x <- seq(from=-1, to=8, by=0.01)
> null.plot <- dt(x, df=100)
> plot(null.plot ~ x, type="n", ylab="probability density", xlab="estimated t-statistic",
+      ylim=c(0, 0.55))
> 
> x.shade <- seq(from=crit, to=8, by=0.01)
> y.shade <- dt(x.shade, df=100, ncp=3)
> polygon(c(min(x.shade),x.shade), c(min(y.shade), y.shade), col=gray(0.8, alpha=1), border=NA)
> 
> y.shade <- dt(x.shade, df=100)
> polygon(c(min(x.shade),x.shade), c(min(y.shade), y.shade), col=gray(0.5, alpha=1), border=NA)
> 
> null.plot <- dt(x, df=100)
> lines(null.plot ~ x)
> 
> alt.plot <- dt(x, df=100, ncp=3)
> lines(alt.plot ~ x, lwd=2)
> 
> abline(v = crit, lty=2)
> 
> arrows(5.5, 0.3, 3.75, 0.15, length=0.1)
> power.value <- round( pt(crit, df=100, ncp=3, lower.tail = F) + pt(-crit, df=100, ncp=3, lower.tail = T), 4)
> text(6, 0.325, labels=(paste( c("power = ",  power.value), collapse = " "  ) ) )
> 
> arrows(5.25, 0.15, 3, 0.025, length=0.1)
> size.value <- round( pt(crit, df=100, lower.tail = F) + pt(-crit, df=100, lower.tail = T), 4)
> text(6.5, 0.175, labels=(paste( c("size = ", size.value), collapse = " "  ) ) )
> 
> legend("topright", lty=c(1,1,2), lwd=c(1,2,1), pt.cex=1, cex=0.8,
+        legend=c("t density under null", "t density, ncp = 3", "critical-t, alpha = 0.005"))
> 
> dev.off()
null device 
          1 
> 
> 
> 
> 
> 
> 
> 
> rm(list=ls())
> set.seed(123456)
> 
> png(filename = "power-comparison.png", width=6, height=6.5, units="in", res=400)
> 
> size <- seq(from=0.5, to=0.9975, by=.0001)
> crit <- qt(size, df=100)
> power <- pt(crit, df=100, ncp=3, lower.tail = F) + pt(-1*crit, df=100, ncp=3, lower.tail = T)
> plot(power ~ size, type="l", ylim=c(0,1))
> 
> size <- seq(from=0.5, to=0.9975, by=.0001)
> crit <- qt(size, df=100)
> power <- pt(crit, df=100, ncp=2, lower.tail = F) + pt(-1*crit, df=100, ncp=2, lower.tail = T)
> lines(power ~ size, type="l", lty=2)
> 
> size <- seq(from=0.5, to=0.9975, by=.0001)
> crit <- qt(size, df=100)
> power <- pt(crit, df=100, ncp=1, lower.tail = F) + pt(-1*crit, df=100, ncp=1, lower.tail = T)
> lines(power ~ size, type="l", lwd=2, col="gray")
> 
> legend("bottomleft", lty=c(1,2,1), lwd=c(1,1,2), col=c("black", "black", "gray"),
+        legend=c("beta = 3", "beta = 2", "beta = 1"))
> 
> dev.off()
null device 
          1 
> 
> 
> 
> png(filename = "multiple-test-power.png", width=6, height=6.5, units="in", res=400)
> 
> coef <- seq(from=0.01, to=6, by=0.01)
> crit <- qt(0.9975, df=100)
> power <- pt(crit, df=100, ncp=coef, lower.tail = F) + pt(-1*crit, df=100, ncp=coef, lower.tail = T)
> crit.2 <- qt( 1-(0.005^(1/2))/2, df=100)
> power.2 <- pt(crit.2, df=100, ncp=coef, lower.tail = F) + pt(-1*crit.2, df=100, ncp=coef, lower.tail = T)
> power.2 <- power.2^2
> crit.3 <- qt( 1-(0.005^(1/3))/2, df=100)
> power.3 <- pt(crit.3, df=100, ncp=coef, lower.tail = F) + pt(-1*crit.3, df=100, ncp=coef, lower.tail = T)
> power.3 <- power.3^3
> 
> plot(power.3 ~ coef, type="l", lwd=2, col="gray", ylim=c(0,1),
+      ylab = "power", xlab = expression(paste("true relationship ", beta)))
> lines(power ~ coef)
> lines(power.2 ~ coef, lty=2)
> 
> legend("bottomright", lty=c(1,2,1), lwd=c(1,1,2),
+        col=c("black", "black", "gray"), 
+        legend = c("one test", "two tests", "three tests"))
> 
> dev.off()
null device 
          1 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> # make a plot of bias associated with lowered statistical significance
> # for a specific relationship, beta = 3, sd_hat-beta = 1, df = 100
> # some code used from: https://www.r-bloggers.com/setting-graph-margins-in-r-using-the-par-function-and-lots-of-cow-milk/
> 
> 
> rm(list=ls())
> set.seed(123456)
> 
> png(filename = "pub-bias-single.png", width=6, height=6.5, units="in", res=400)
> 
> alpha.v <- seq(from=0.005, to=0.1, by=0.001)
> bias.v <- rep(NA, length(alpha.v))
> crit.v <- qt( 1-(alpha.v/2), df=100 )
> 
> samples <- rt(1e7, ncp=3, df=100)
> 
> pb <- txtProgressBar(min=1, max=length(crit.v), initial=0, style=3)
> for(i in 1:length(crit.v)){
+   
+   setTxtProgressBar(pb, i)
+   bias.v[i] <- mean( (samples[ which(samples > crit.v[i] ) ]  / 3) - 1 )
+   
+ }
  |                                                                              |                                                                      |   0%  |                                                                              |=                                                                     |   1%  |                                                                              |=                                                                     |   2%  |                                                                              |==                                                                    |   3%  |                                                                              |===                                                                   |   4%  |                                                                              |====                                                                  |   5%  |                                                                              |====                                                                  |   6%  |                                                                              |=====                                                                 |   7%  |                                                                              |======                                                                |   8%  |                                                                              |=======                                                               |   9%  |                                                                              |=======                                                               |  11%  |                                                                              |========                                                              |  12%  |                                                                              |=========                                                             |  13%  |                                                                              |==========                                                            |  14%  |                                                                              |==========                                                            |  15%  |                                                                              |===========                                                           |  16%  |                                                                              |============                                                          |  17%  |                                                                              |=============                                                         |  18%  |                                                                              |=============                                                         |  19%  |                                                                              |==============                                                        |  20%  |                                                                              |===============                                                       |  21%  |                                                                              |===============                                                       |  22%  |                                                                              |================                                                      |  23%  |                                                                              |=================                                                     |  24%  |                                                                              |==================                                                    |  25%  |                                                                              |==================                                                    |  26%  |                                                                              |===================                                                   |  27%  |                                                                              |====================                                                  |  28%  |                                                                              |=====================                                                 |  29%  |                                                                              |=====================                                                 |  31%  |                                                                              |======================                                                |  32%  |                                                                              |=======================                                               |  33%  |                                                                              |========================                                              |  34%  |                                                                              |========================                                              |  35%  |                                                                              |=========================                                             |  36%  |                                                                              |==========================                                            |  37%  |                                                                              |===========================                                           |  38%  |                                                                              |===========================                                           |  39%  |                                                                              |============================                                          |  40%  |                                                                              |=============================                                         |  41%  |                                                                              |=============================                                         |  42%  |                                                                              |==============================                                        |  43%  |                                                                              |===============================                                       |  44%  |                                                                              |================================                                      |  45%  |                                                                              |================================                                      |  46%  |                                                                              |=================================                                     |  47%  |                                                                              |==================================                                    |  48%  |                                                                              |===================================                                   |  49%  |                                                                              |===================================                                   |  51%  |                                                                              |====================================                                  |  52%  |                                                                              |=====================================                                 |  53%  |                                                                              |======================================                                |  54%  |                                                                              |======================================                                |  55%  |                                                                              |=======================================                               |  56%  |                                                                              |========================================                              |  57%  |                                                                              |=========================================                             |  58%  |                                                                              |=========================================                             |  59%  |                                                                              |==========================================                            |  60%  |                                                                              |===========================================                           |  61%  |                                                                              |===========================================                           |  62%  |                                                                              |============================================                          |  63%  |                                                                              |=============================================                         |  64%  |                                                                              |==============================================                        |  65%  |                                                                              |==============================================                        |  66%  |                                                                              |===============================================                       |  67%  |                                                                              |================================================                      |  68%  |                                                                              |=================================================                     |  69%  |                                                                              |=================================================                     |  71%  |                                                                              |==================================================                    |  72%  |                                                                              |===================================================                   |  73%  |                                                                              |====================================================                  |  74%  |                                                                              |====================================================                  |  75%  |                                                                              |=====================================================                 |  76%  |                                                                              |======================================================                |  77%  |                                                                              |=======================================================               |  78%  |                                                                              |=======================================================               |  79%  |                                                                              |========================================================              |  80%  |                                                                              |=========================================================             |  81%  |                                                                              |=========================================================             |  82%  |                                                                              |==========================================================            |  83%  |                                                                              |===========================================================           |  84%  |                                                                              |============================================================          |  85%  |                                                                              |============================================================          |  86%  |                                                                              |=============================================================         |  87%  |                                                                              |==============================================================        |  88%  |                                                                              |===============================================================       |  89%  |                                                                              |===============================================================       |  91%  |                                                                              |================================================================      |  92%  |                                                                              |=================================================================     |  93%  |                                                                              |==================================================================    |  94%  |                                                                              |==================================================================    |  95%  |                                                                              |===================================================================   |  96%  |                                                                              |====================================================================  |  97%  |                                                                              |===================================================================== |  98%  |                                                                              |===================================================================== |  99%  |                                                                              |======================================================================| 100%> close(pb)

> 
> bias.v <- 100*bias.v
> 
> plot(bias.v ~ alpha.v, type="l", xlim=c(0, 0.1), ylim=c(0, 30),
+      ylab = "estimated publication bias (as % of coefficient)",
+      xlab = "alpha for statistical significance test")
> 
> dev.off()
null device 
          1 
> 
> # generate a plot for an overall prior distribution including a spike
> # at beta = 0
> 
> rm(list=ls())
> set.seed(123456)
> 
> png(filename = "pub-bias-dist.png", width=6, height=6.5, units="in", res=400)
> 
> gen.norm.samp <- function(se.b, df.b, pr.false, n.draws, bound, seed){
+   
+   if(seed!="off"){
+     set.seed(seed)
+   }
+   b.true <- rnorm(n.draws, mean=0, sd=bound)
+   b.false <- rep(0, n.draws)
+   choice <- ifelse(runif(n.draws)<(1-pr.false), 1, 0)
+   b.test <- ifelse(choice==1, b.true, b.false)
+   normalizer <- mean(abs(b.test))
+   
+   b.test.est <- b.test + se.b*rt(n.draws, df=df.b)
+   bias <- b.test.est - b.test
+   
+   mult <- sign(b.test)
+   mult <- ifelse(mult==0, sign(bias), mult)
+   bias.out <- (mult*bias)/normalizer
+   
+   b.t <- abs(b.test.est/se.b)
+   
+   list.out <- list(b.true = b.test, b.est = b.test.est, b.bias = bias.out, b.t = b.t)
+   return(list.out)
+   
+ }
> 
> 
> alpha.v <- seq(from=0.005, to=0.1, by=0.001)
> bias.v <- rep(NA, length(alpha.v))
> crit.v <- qt( 1-(alpha.v/2), df=100 )
> 
> samples.list <- gen.norm.samp(se.b=1, df.b=100, pr.false=0.5, n.draws=1e7, bound=3, seed=123456)
> 
> pb <- txtProgressBar(min=1, max=length(crit.v), initial=0, style=3)
> for(i in 1:length(crit.v)){
+   
+   setTxtProgressBar(pb, i)
+   sig <- which(samples.list$b.t > crit.v[i])
+   bias.v[i] <- mean(samples.list$b.bias[sig])
+   
+ }
  |                                                                              |                                                                      |   0%  |                                                                              |=                                                                     |   1%  |                                                                              |=                                                                     |   2%  |                                                                              |==                                                                    |   3%  |                                                                              |===                                                                   |   4%  |                                                                              |====                                                                  |   5%  |                                                                              |====                                                                  |   6%  |                                                                              |=====                                                                 |   7%  |                                                                              |======                                                                |   8%  |                                                                              |=======                                                               |   9%  |                                                                              |=======                                                               |  11%  |                                                                              |========                                                              |  12%  |                                                                              |=========                                                             |  13%  |                                                                              |==========                                                            |  14%  |                                                                              |==========                                                            |  15%  |                                                                              |===========                                                           |  16%  |                                                                              |============                                                          |  17%  |                                                                              |=============                                                         |  18%  |                                                                              |=============                                                         |  19%  |                                                                              |==============                                                        |  20%  |                                                                              |===============                                                       |  21%  |                                                                              |===============                                                       |  22%  |                                                                              |================                                                      |  23%  |                                                                              |=================                                                     |  24%  |                                                                              |==================                                                    |  25%  |                                                                              |==================                                                    |  26%  |                                                                              |===================                                                   |  27%  |                                                                              |====================                                                  |  28%  |                                                                              |=====================                                                 |  29%  |                                                                              |=====================                                                 |  31%  |                                                                              |======================                                                |  32%  |                                                                              |=======================                                               |  33%  |                                                                              |========================                                              |  34%  |                                                                              |========================                                              |  35%  |                                                                              |=========================                                             |  36%  |                                                                              |==========================                                            |  37%  |                                                                              |===========================                                           |  38%  |                                                                              |===========================                                           |  39%  |                                                                              |============================                                          |  40%  |                                                                              |=============================                                         |  41%  |                                                                              |=============================                                         |  42%  |                                                                              |==============================                                        |  43%  |                                                                              |===============================                                       |  44%  |                                                                              |================================                                      |  45%  |                                                                              |================================                                      |  46%  |                                                                              |=================================                                     |  47%  |                                                                              |==================================                                    |  48%  |                                                                              |===================================                                   |  49%  |                                                                              |===================================                                   |  51%  |                                                                              |====================================                                  |  52%  |                                                                              |=====================================                                 |  53%  |                                                                              |======================================                                |  54%  |                                                                              |======================================                                |  55%  |                                                                              |=======================================                               |  56%  |                                                                              |========================================                              |  57%  |                                                                              |=========================================                             |  58%  |                                                                              |=========================================                             |  59%  |                                                                              |==========================================                            |  60%  |                                                                              |===========================================                           |  61%  |                                                                              |===========================================                           |  62%  |                                                                              |============================================                          |  63%  |                                                                              |=============================================                         |  64%  |                                                                              |==============================================                        |  65%  |                                                                              |==============================================                        |  66%  |                                                                              |===============================================                       |  67%  |                                                                              |================================================                      |  68%  |                                                                              |=================================================                     |  69%  |                                                                              |=================================================                     |  71%  |                                                                              |==================================================                    |  72%  |                                                                              |===================================================                   |  73%  |                                                                              |====================================================                  |  74%  |                                                                              |====================================================                  |  75%  |                                                                              |=====================================================                 |  76%  |                                                                              |======================================================                |  77%  |                                                                              |=======================================================               |  78%  |                                                                              |=======================================================               |  79%  |                                                                              |========================================================              |  80%  |                                                                              |=========================================================             |  81%  |                                                                              |=========================================================             |  82%  |                                                                              |==========================================================            |  83%  |                                                                              |===========================================================           |  84%  |                                                                              |============================================================          |  85%  |                                                                              |============================================================          |  86%  |                                                                              |=============================================================         |  87%  |                                                                              |==============================================================        |  88%  |                                                                              |===============================================================       |  89%  |                                                                              |===============================================================       |  91%  |                                                                              |================================================================      |  92%  |                                                                              |=================================================================     |  93%  |                                                                              |==================================================================    |  94%  |                                                                              |==================================================================    |  95%  |                                                                              |===================================================================   |  96%  |                                                                              |====================================================================  |  97%  |                                                                              |===================================================================== |  98%  |                                                                              |===================================================================== |  99%  |                                                                              |======================================================================| 100%> close(pb)

> 
> par(mar=c(5,6,4,2)+0.1,mgp=c(4,1,0))
> 
> bias.v <- 100*bias.v
> 
> # add line to the plot
> plot(bias.v ~ alpha.v, type="l", xlim=c(0, 0.1), ylim=c(40, 50),
+      ylab = c("estimated publication bias","(as % of mean relationship in prior)"),
+      xlab = " ")
> mtext(text = "alpha for statistical significance test", side=1, line=3)
> 
> dev.off()
null device 
          1 
> 
> 
> proc.time()
   user  system elapsed 
  38.89   11.10   51.32 
